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A comprehensive survey of event-by-event fluctuations of charged hadron multiplicity in rela- 
tivistic heavy ions is presented. The survey covers Au+Au collisions at -^snn ~ 62.4 and 200 
GeV, and Cu-I-Cu collisions at .^snn = 22.5, 62.4, and 200 GeV. Fluctuations are measured as a 
function of collision centrality, transverse momentum range, and charge sign. After correcting for 
non-dynamical fluctuations due to fluctuations in the collision geometry within a centrality bin, the 
remaining dynamical fluctuations expressed as the variance normalized by the mean tend to decrease 
with increasing centrality. The dynamical fluctuations are consistent with or below the expectation 
from a superposition of participant nucleon-nucleon collisions based upon p-|-p data, indicating that 
this dataset does not exhibit evidence of critical behavior in terms of the compressibility of the 
system. A comparison of the data with a model where hadrons are independently emitted from a 
number of hadron clusters suggests that the mean number of hadrons per cluster is small in heavy 
ion collisions. 

PACS numbers: 25.75.Gz, 25.75.Nq, 21.65.Qr, 25.75.Ag 



I. INTRODUCTION 

Recent work with lattice gauge theory simulations 
has attempted to map out the phase diagram of Quan- 
tum Chromodynamics (QCD) in temperature and baryo- 
chemical potential (/ie) using finite values of the up and 
down quark masses. The results of these studies indi- 
cate that the QCD phase diagram may contain a first- 
order transition line between the hadron gas phase and 
the strongly-coupled Quark-Gluon Plasma (sQGP) phase 
that terminates at a critical point This property 
is analogous to that observed in the phase diagram for 
many common liquids and other substances, including 
water. However, different model predictions and lattice 
calculations yield widely varying estimates of the loca- 
tion of the critical point on the QCD phase diagram 
0. Direct experimental observation of critical phenom- 
ena in heavy ion collisions would confirm the existence of 
the critical point, narrow down its location on the QCD 
phase diagram, and provide an important constraint for 
the QCD models. 

The estimated value of energy densities achieved in 
heavy ion collisions at the Brookhaven National Labo- 
ratory's Relativistic Heavy Ion Collider (RHIC) exceeds 



'Deceased 

tPHENIX Spokesperson: j acak@skipper .physics . suny sb . edu | 



the threshold for a phase transition from normal hadronic 
matter to partonic matter. Recent experimental evidence 
indicates that properties of the matter being produced 
include strong collective flow and large opacity to scat- 
tered quarks and gluons - the matter appears to behave 
much like a perfect fluid Q . While measurements suggest 
the produced matter has properties that differ from nor- 
mal nuclear matter, unambiguous evidence of the nature 
and location of any phase transition from normal nuclear 
matter has been elusive thus far. Described here is a 
search for direct evidence of a phase transition by mea- 
suring the fluctuations of the event-by-event multiplici- 
ties of produced charge particles in a variety of collision 
systems. 

In order to illustrate how the measurement of charged 
particle multiplicity fluctuations can be sensitive to the 
presence of a phase transition, the isothermal compress- 
ibility of the system can be considered Q . The isother- 
mal compressibility is defined as follows: 

kT^-l/V{SV/dP)^, (1) 

where V is the volume, T is the temperature, and P is 
the pressure of the system. In order to relate the com- 
pressibility to measurements of multiplicity fluctuations, 
we assume that relativistic nucleus-nucleus collisions can 
be described as a thermal system in the Grand Canon- 
ical Ensemble (GCE) Q. The GCE can be applied to 
the case of measurements near mid-rapidity since energy 
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and conserved quantum numbers in this region can be ex- 
changed with the rest of the system, that serves as a heat 
bath [1] . Detailed studies of multipUcity fluctuations in 
the Canonical and Microcanonical Ensembles with the 
application of conservation laws can be found elsewhere 
P, H]- In the GCE, the isothermal compressibility is di- 
rectly related to the variance of the particle multiplicity 
as follows: 



{{N - {N)f) =var{N) = 



V 



(2) 



where N is the particle multiplicity, (N) = /i^r is 
the mean multiplicity, and fee is Boltzmann's constant 
0. Here, multiplicity fluctuation measurements are pre- 
sented in terms of the scaled variance, u-^: 



var{N) UN , 



(3) 



In a continuous, or second-order, phase transition, the 
compressibility diverges to an infinite value at the critical 
point. Near the critical point, this divergence is described 
by a power law in the variable e — {T — Tq) /Tq, where 
Tc is the critical temperature. Hence, the relationship 
between multiplicity fluctuations and the compressibility 
can be exploited to search for a clear signature of critical 
behavior by looking for the expected power law scaling 
of the compressibility: 



«T oc (— ) ^ cx e 

In 



(4) 



where 7 is the critical exponent for isothermal compress- 
ibility [9|]. If the QCD phase diagram contains a critical 
point, systems with a low value of baryo-chemical poten- 
tial (/is) could pass through the cross-over region and un- 
dergoa continuous phase transition [2| . Recent estimates 
[lol of the behavior of the quark number susceptibil- 
ity, Xq, which is proportional to the value of the isother- 
mal compressibility of the system, predict that its value 
will increase by at least an order of magnitude close to the 
QCD critical point. Given that the scaled variance is pro- 
portional to fcx, measurements of charged particle mul- 
tiplicity are expected to be a sensitive probe for critical 
behavior. In addition, within a scenario where droplets 
of Quark-Gluon Plasma are formed during a first-order 
phase transition, the scaled variance of the multiplicity 
could increase by a factor of 10-100 [121] . 

Experimentally, a search for critical behavior is facili- 
tated by the rich and varied dataset provided by RHIC. 
It is expected that the trajectory of the colliding sys- 
tem in the QCD phase diagram can be modified by vary- 
ing the colliding energy If the system approaches 
close enough to the critical line for a long enough time 
period, then critical phenomena could be readily appar- 
ent through the measurement of multiplicity fluctuations 
[1^. It may also be possible to determine the critical 
exponents of the system. Nature tends to group ma- 
terials into universality classes whereby all materials in 



the same universality class share identical values for their 
set of critical exponents. Although beyond the scope of 
this analysis, observation of critical behavior in heavy ion 
collisions and the subsequent measurement of the critical 
exponents could determine the universality class in which 
QCD is grouped, providing essential constraints for the 
models. 



Charged particle multiplicity fluctuations have been 
measured in elementary collisions over a large range of 
collisions energies fll, [ll, [H, [13, [H, [H, [13 ■ Initial mea- 
surements of multiplicity fluctuations in minimum-bias 
0+Cu colhsions at v^Snn=4.86 GeV were made by BNL 
Experiment E802 [2l|, minimum-bias 0-1- Au collisions 
at y5ii7=17.3 GeV by CERN Experiment WA80 [H 
and minimum-bias S-f S, 0-|-Au, and S-|-Au collisions at 
7^=17.3 GeV by CERN Experiment NA35 [H. Re- 
cently, larger datasets have enabled the measurement of 
the centrality-dependence of multiplicity fluctuations in 
Pb-fPb collisions at ^3^^7=17.3 GeV by CERN Exper- 
iment WA98 HI and in Pb-hPb, C-fC, and Si-hSi col- 
lisions at ysii7=17.3 GeV by CERN Experiment NA49 
[25]. The PHENIX Experiment at RHIC has performed 
an analysis of density correlations in longitudinal space 
with a differential analysis of charged particle multiplic- 
ity fluctuations in 200 GeV Au-|-Au collisions over the 
entire transverse momentum range [2^. Thus far, the 
fluctuation measurements in heavy ion collisions do not 
indicate significant signs of a phase transition. However, 
the full range of collision energies and species accessible 
by RHIC are yet to be explored. 



Presented here is a comprehensive survey of multi- 
plicity fluctuations of charged hadrons measured by the 
PHENIX Experiment at RHIC. The survey will cover 
the following coUision systems: ^snn=200 GeV Au-|-Au, 
62.4 GeV Au-hAu, 200 GeV Cu-f Cu, 62.4 GeV Cu-I-Cu, 
and 22.5 GeV Cu-|-Cu with comparisons to -^5=200 GeV 
p-|-p collisions, which serve as a baseline measurement. 
The Au-I-Au data were taken during RHIC Run-4 (2004), 
the Cu-f-Cu data were taken during RHIC Run-5 (2005), 
and the p+p data were taken during RHIC Run-3 (2003). 
Multiplicity fluctuations for each collision system with 
the exception of p+p will also be presented as a function 
of centrality to help select the system volume. Multiplic- 
ity fluctuations will also be presented as a function of 
transverse momentum range, and charge sign. 



This paper is organized as follows: Sec. II will discuss 
the experimental apparatus and details; Sec. Ill will dis- 
cuss the methods applied for the measurement of mul- 
tiplicity fluctuations and the removal of non-dynamical 
fluctuations due to fluctuations of the collision geometry 
within a centrality bin; Sec. IV will present the results 
and compare them to other models. Sec. V will present 
a discussion and summary of the results. 
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II. EXPERIMENTAL SETUP 

The PHENIX detector consists of two central spec- 
trometer arms designed for charged particle tracking, 
designated east and west, and two muon spectrometers 
designed for muon tracking and identification, designated 
north and south. The muon spectrometers are not used 
in this analysis. A comprehensive description of the 
PHENIX detector is documented elsewhere [l^l- The 
analysis described here utilizes the central spectrometer 
arms, which consist of a set of tracking detectors [28| . 
particle identification detectors [2§|, and an electromag- 
netic calorimeter |3C|. The central spectrometer arms 
cover a rapidity range of |r/| < 0.35 and each arm sub- 
tends 90 degrees in azimuth. A detailed description of the 
algorithms and performance of the central arm track re- 
consti'uction and momentum reconstruction can be found 
in (3l|. 

There are two detectors that are used for triggering, 
centrality determination, and event vertex determina- 
tion. The Beam-Beam Counters (BBCs) consist of 64 
individual quartz Cherenkov counters that cover the full 
azimuthal angle in the pseudorapidity range 3.0 < I77I < 
3.9. The Zero Degree Calorimeters (ZDCs) cover the 
pseudorapidity range I77I > 6 and measure the energy of 
spectator neutrons with an energy resolution of approx- 
imately 20%. More details about these detectors can be 
found in [3^ . The collision vertex position is determined 
using timing information from the BBCs with an r.m.s. 
resolution for central Au-|-Au events of 6 mm along the 
beam axis. The collision vertex is required to be recon- 
structed within ±30 cm from the center of the spectrom- 
eter. The BBCs also provide a minimum biased (MB) 
event trigger. 

Due to the large dynamic range in -^/snn covered by 
this analysis, it is necessary to implement algorithms that 
are dependent on the collision energy for the determina- 
tion of the centrality of each event. In Au+Au collisions 
at -y/SNN=200 GeV, the centrality of the collision is deter- 
mined by using correlations of the total energy deposited 
in the ZDCs with the total charge deposited in the BBCs 
as described in However, in 200 GeV Cu-hCu, 62.4 
GeV Cu-I-Cu, and 62.4 GeV Au-I-Au colhsions, the re- 
solving power of the ZDCs is insufficient to significantly 
contribute to the centrality definition. Therefore, only 
the total charge deposited in the BBCs is used to deter- 
mine centrality in these collision systems, as described in 
[3^. Using the 200 GeV Au-I-Au data, it has been veri- 
fied that application of the BBC-ZDC correlation for the 
centrality definition as opposed to the BBC-only defini- 
tion shows no significant differences in the values of the 
charged hadron fiuctuation quantities presented here as 
a function of centrality. 

The location of the BBCs are fixed for every collision 
energy. At the lowest collision energy (^snn=22.5 GeV), 
it becomes kinematically possible for spectator nucleons 
to fall within the acceptance of the BBC. This results 
in a BBC response in its total charge sum that is no 



longer linear with the number of participating nucleons 
(A^part)- In this case, it becomes necessary to define the 
centrality using the total charged particle multiplicity in 
Pad Chamber 1 (PCI) PCI is chosen due to its 

fine segmentation, high tracking efficiency, and relative 
proximity to the event vertex. Details on this procedure 
are also described in ^33\. For all collision species and 
energies, the distribution of the number of participants 
was determined using a Monte Carlo simulation based 
upon the Glauber model [H, [s^l • 

The number of minimum bias events analyzed for each 
dataset are 25.6 million events for 200 GeV Au-|-Au, 24.9 
million events for 62.4 GeV Au-|-Au, 15.0 million events 
for 200 GeV Cu+Cu, 12.2 million events for 62.4 GeV 
Cu-I-Cu, 5.5 million events for 22.5 GeV Cu+Cu, and 2.7 
million events for 200 GeV p-l-p. Only a fraction of the 
complete 200 GeV Au-|-Au, Cu-f Cu, and p-fp datasets 
are analyzed, but this fraction is more than sufficient for 
this analysis. 

The charged particle multiplicity is determined on an 
event-by-event basis by counting the number of unam- 
biguous reconstructed tracks in the Drift Chamber origi- 
nating from the collision vertex that have corresponding 
hits in Pad Chamber 1 and Pad Chamber 3. Track se- 
lection includes cuts on reconstructed tracks in the Drift 
Chamber to reduce double-counted ghost tracks to a neg- 
ligible level. In order to minimize background originat- 
ing from the magnets, reconstructed tracks are required 
to lie within ±75 cm from the center of the Drift Cham- 
ber along the beam axis. This requirement reduces the 
pseudorapidity range of reconstructed tracks to \ri\ <0.26. 
The Ring Imaging Cherenkov detector (RICH) is utilized 
to reduce background from electrons resulting from pho- 
ton conversions. 

Although the central arm spectrometer covers a to- 
tal azimuthal range of tt radians, detector and tracking 
inefficiencies reduce the effective average azimuthal ac- 
tive area to 2.1 radians for the 200 Gev Au+Au and 200 
GeV p+p datasets, and 2.0 radians for the other datasets. 
Fluctuation quantities are quoted for these acceptances 
separately for each dataset. The differences in accep- 
tance between datasets, which are due to variations in 
the detector over the three year period in which the data 
was collected, result in less than a 1% variation in the 
fluctuation quantities quoted here. 

III. DATA ANALYSIS 

Multiplicity fluctuations of charged particles, desig- 
nated Wch, can be generally defined [35| as follows: 

^ch — 7T^-^ — , [O) 

where Nch is the charged particle multiplicity. Simply 
stated, the fluctuations can be quoted as the variance 
of the multiplicity (cr^jj) normalized by the mean {/ich = 
( A^'ch) ) • This is also referred to as the scaled variance [1^ . 
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If the multiplicity distribution is Poissonian, the scaled 
variance is 1.0. 

It has been well established that charged particle mul- 
tiplicity distributions in elementary nucleon-nucleon col- 
lisions can be described by the Negative Binomial Distri- 
bution (NBD) [U in, [3 . The NBD also well describes 
multiplicity distributions in heavy ion collisions [U, [22| . 
The Negative Binomial Distribution of an integer n is 
defined as follows: 



P{n) 



r(n + fcNBo) 



(Mch/fcNBD)" 



r(n + l)r(A:NBD) (1 + Mch/fcNBD)"+'=^ 



(6) 



where P(n) is normalized to 1.0 over the range < n < 
oo, /ich = (Nch) = (n), and /cnbd is an additional pa- 
rameter. The NBD reduces to a Poisson distribution in 
the limit /cnbd ^ oo. The NBD variance and mean is 
related to /cnbd as follows: 



^ch 



1 



1 



k 



NBD 



Hence, the scaled variance is given by 



1 



Mch 



BD 



(7) 



(8) 



A useful property of the Negative Binomial Distribu- 
tion concerns its behavior when a population that fol- 
lows the NBD is subdivided randomly by repeated inde- 
pendent trials with a constant probability onto smaller 
subsets. This results in a binomial decomposition of the 
original population into subsets that also follow the NBD 
with the same value of /cnbd [21| . This property can be 
applied to estimate the behavior of multiplicity fluctu- 
ations as a function of acceptance, assuming that there 
are no significant correlations present over the acceptance 
range being examined. Starting with an original NBD 
sample with mean /Zch and scaled variance Wch, a sam- 
ple in a fractional acceptance with mean /iacc is also de- 
scribed by an NBD distribution. An acceptance fraction 
can be defined as /acc = /-tacc/z^ch- The scaled variance 
of the subsample from Equation ^ is thus 

Waco = 1 + (Aiacc/feNBo) = 1 + (./accAich/fcNBD ) • (9) 

Since /cnbd is identical for the two samples, /ich/^NBD = 
Wch — 1 can be substituted, yielding the following relation 
between the scaled variances of the original and fractional 
acceptance samples: 



1 + /acc('^ch - 1) 



(10) 



Thus, the measured scaled variance will decrease as the 
acceptance is decreased while /cnbd remains constant, if 
there are no additional correlations present over the given 
acceptance range. 

Figures [T][2] show the uncorrected, or raw, multiplicity 
distributions in the pr range 0.2 < pT < 2.0 GeV for all 
centralities from each collision system overlayed with fits 
to Negative Binomial Distributions (dashed lines). For 



presentation purposes, the data have been normalized 
on the horizontal axis by the mean of the distribution 
and scaled on the vertical axis by the successive amounts 
stated in the legend. The NBD fits describe the data 
distributions very well for all collision systems, centrali- 
ties, and pt ranges. Hence, the mean and variance of the 
multiplicity distributions presented here are all extracted 
from NBD fits. The results of each fit for 0.2 < pr < 2.0 
GeV are compiled in Table HI The mean and standard 
deviation of each fit for 0.2 < px < 2.0 GeV are plotted 
in Fig. O 

Each dataset was taken over spans of several days to 
several weeks, all spanning three separate RHIC run- 
ning periods. During these periods, changes in the total 
acceptance and efficiency of the central arm spectrom- 
eters cause the fluctuation measurements to vary, thus 
introducing an additional systematic error to the results. 
This systematic error was minimized by requiring that 
the dataset is stable in quantities that are sensitive to 
detector variations, including the mean charged parti- 
cle multiplicity, mean collision vertex position, and mean 
centrality. A time-dependent systematic error is applied 
independently to each point by calculating the standard 
deviation of the scaled variance calculated from subsets 
of the entire dataset, with each subset containing about 
1 million events. These systematic errors are applied to 
all subsequent results. 

The tracking efficiency of the PHENIX central arm 
spectrometer is dependent on centrality, especially in the 
most central 200 GeV Au-|-Au collisions [30|. With the 
assumption that tracking inefhciencies randomly effect 
the multiplicity distribution on an event-by-event basis, 
the effect of inefficiencies on the scaled variance can be 
estimated using Equation pH)) where /acc is replaced by 
the inverse of the tracking efficiency, l//cff. Tracking effi- 
ciency effects the value of the scaled variance by 1.5% at 
the most. The scaled variance has been corrected for 
tracking inefficiency as a function of centrality for all 
species. The uncertainty of the tracking efficiency es- 
timate is typically 2% and has been propagated into the 
systematic error estimate on a point-by-point basis. 

Due to the non-zero width of the centrality bin selec- 
tion from the data, each centrality bin necessarily se- 
lects a range of impact parameters. This introduces a 
non-dynamical fluctuation component to the measured 
multiplicity fluctuations due to the resulting fluctuations 
in the geometry of the collisions [H, [37|]. Therefore, it is 
necessary to estimate the magnitude of the geometry fluc- 
tuation component so that only the interesting dynami- 
cal fluctuations remain. The most practical method for 
estimating the geometry fluctuation component is with 
a model of heavy ion collisions. The URQMD [3^ and 
HSD 40, 41] models have previously been applied for this 
purpose. Here, the HIJING event generator [s^l is cho- 
sen for this estimate because it well reproduces the mean 
multiplicity in heavy ion collisions |33| as measured by 
the PHENIX detector. HIJING includes multiple mini- 
jet production based upon QCD-inspired models, soft ex- 
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citation, nuclear shadowing of parton distribution func- 
tions, and the interaction of jets in dense nuclear matter. 
The estimate is performed individually for each central- 
ity bin, collision system, and pt range using the following 
procedure. First, HIJING is run with an impact param- 
eter distribution that is sampled from a Gaussian distri- 
bution with a mean and standard deviation that, for a 
given centrality bin, reproduces the distributions of the 
charge deposited in the BBC and the energy deposited in 
the ZDC (for 200 GeV Au-hAu). Second, HIJING is run 
at a fixed impact parameter with a value identical to the 
mean of the Gaussian distribution in the first run. For 
each centrality bin, 12,000 HIJING events are processed 
for each impact parameter selection. The scaled variance 
for each impact parameter selection, cjcauss and CLifixcd, is 
extracted and the measured scaled variance is corrected 
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FIG. 1: The uncorrected multiplicity distributions of charged 
hadrons with 0.2 < pr < 2.0 GeV/c for 200 (upper) and 62.4 
(lower) GeV Au-|-Au collisions. The dashed lines are fits to 
the Negative Binomial Distribution. The data are normalized 
to the mean and scaled by the amounts in the legend. 



FIG. 2: The uncorrected multiplicity distributions of charged 
hadrons with 0.2 < pr < 2.0 GeV/c for 200 (upper), 62.4 
(middle), and 22.5 (lower) GeV Cu-fCu and 200 GeV p+p 
(lower) collisions. The dashed lines are fits to the Negative 
Binomial Distribution. The data are normalized to the mean 
and scaled by the amounts in the legend. 
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TABLE I: Tabulation of the charged hadron multiphcity data and corrections for 0.2 < pr < 2.0 GeV/c. The errors quoted for 
/ich and (Jch represent their time-dependent systematic error. The errors quoted for a;ch,dyn and l/feNBD,dyn represent their total 
systematic error. For each dataset the first three columns give the species, collision energy, and geometric correction factor, 

/gco, respectively. 



Species 


(GeV) 


/geo 


A'part 


A*ch 


raw cr^h 


^ch,dyn 




l/fcNBD,dyn 




xVdof 








351 


61.0 ± 1.1 


75.6 ± 1.9 


1.10 ±0.02 


1.45 


10"°^ ± 2.2 • 




37.1/58 








299 


53.1 ± 1.0 


71.8 ±1.8 


1.15 ±0.02 


2.45 


10""^ ±2.7- 


10"°^ 


38.6/56 








253 


45.8 ± 0.8 


65.2 ± 1.5 


1.17 ±0.02 


3.41 


10""^ ± 2.9 • 




34.0/54 








215 


39.1 ±0.7 


57.8 ± 1.6 


1.19 ±0.03 


4.53 


10"°^ ± 3.6 • 


10-04 


29.1/53 








181 


32.6 ±0.6 


49.7 ± 1.3 


1.21 ±0.03 


5.95 


10"°^ ± 5.1 ■ 




24.5/50 


Au+Au 


200 


0.37 ±0.027 


151 


27.4 ± 0.5 


41.4 ± 1.0 


1.20 ± 0.03 


6.86 


10""^ ± 5.5 • 




20.7/46 








125 


22.3 ± 0.4 


33.8 ±0.9 


1.20 ± 0.03 


8.47 






11.9/41 








102 


17.8 ±0.3 


26.7 ±0.6 


1.19 ±0.02 


1.05 


10"°^ ± 9.0 • 


^Q-04 


16.6/37 








82 


14.2 ±0.3 


20.8 ±0.6 


1.17 ±0.02 


1.20 


10""^ ± 1.0 ■ 


10""^ 


37.8/33 








65 


10.8 ±0.2 


16.0 ± 0.4 


1.18 ±0.02 


1.64 


10""^ ± 1.3 • 




37.8/28 








51 


8.3 ±0.2 


12.1 ±0.3 


1.17 ±0.02 


2.06 


10"°^ ± 2.0 • 




53.8/24 








345 


44.0 ±0.3 


53.6 ±0.5 


1.08 ±0.02 


1.63 


10""^ ± 2.0 ■ 


io-"4 


14.6/54 








296 


37.3 ± 0.2 


48.3 ± 0.3 


1.11 ±0.02 


2.63 


10"°^ ± 2.6 • 


^Q-04 


13.8/53 








250 


31.0 ±0.2 


39.8 ±0.4 


1.10 ±0.02 


3.00 


10"°^ ± 3.0 • 


^Q-04 


14.0/50 








211 


25.4 ±0.2 


33.6 ±0.5 


1.12 ±0.02 


4.21 


lO""-'^ ± 4.4 ■ 


^Q-04 


8.36/44 








177 


20.8 ±0.1 


27.8 ±0.2 


1.12 ±0.02 


5.34 


10""'^ ± 5.5 ■ 


^Q-04 


19.2/40 


Au+Au 


62.4 


0.33 ±0.031 


148 


16.6 ±0.1 


22.8 ± 0.3 


1.13 ±0.02 


7.43 


10""^ ± 7.8 • 


^Q-04 


25.9/37 








123 


13.1 ±0.1 


18.1 ±0.2 


1.13 ±0.02 


9.61 


10"°^ ± 9.7 • 


^Q-04 


34.3/33 








102 


10.4 ±0.1 


14.9 ±0.1 


1.15 ±0.02 


1.38 


10"°^ ± 1.4- 




44.5/28 








82 


7.8 ±0.1 


11.1 ±0.1 


1.14 ±0.02 


1.76 


10""^ ± 1.9 • 


2^Q-03 


50.9/24 








66 


5.9 ± 0.04 


8.3 ±0.1 


1.14 ±0.02 


2.37 


10"°^ ± 3.8 • 


j^Q-03 


45.4/20 








51 


4.1 ± 0.03 


5.8 ± 0.04 


1.13 ±0.02 


3.08 


10"°^ ± 9.1 • 


10""^ 


36.2/17 








104 


19.3 ±0.3 


25.7 ± 0.8 


1.14 ±0.03 


6.93 


10"°^ ± 1.3- 


^Q-03 


24.3/30 








92 


16.0 ±0.2 


21.9 ±0.5 


1.15 ±0.03 


9.26 


10"°^ ± 1.5 • 


^Q-03 


21.7/31 








79 


13.5 ±0.2 


18.8 ±0.4 


1.16 ±0.03 


1.15 


10"°^ ± 2.1 • 


^Q-03 


19.4/29 


Cu+Cu 


200 


0.40 ± 0.047 


67 


11.1 ±0.2 


15.3 ±0.3 


1.15 ±0.03 


1.36 


10""^ ± 2.0 ■ 


■^Q-03 


29.9/26 








57 


9.2 ± 0.1 


13.0 ± 0.3 


1.17 ± 0.03 


1.75 




^Q-03 


26.0/25 








48 


7.5 ±0.1 


10.5 ±0.2 


1.16 ±0.03 


2.14 


10-°2 ± 3.6 • 


10-03 


30.6/22 








40 


6.2 ±0.1 


8.7 ±0.2 


1.17 ±0.03 


2.69 


10""^ ±4.8- 


lQ-03 


28.6/20 








33 


4.9 ± 0.06 


6.8 ±0.1 


1.16 ±0.03 


3.12 


10""^ ± 8.5 • 


10-03 


45.7/18 








104 


12.6 ±0.1 


16.7 ±0.2 


1.10 ±0.03 


8.16 


10"°^ ± 1.7- 


10-03 


40.6/31 








92 


11.0 ±0.1 


16.2 ±0.1 


1.15 ±0.04 


1.35 


10"°^ ±2.7- 


10-03 


64.2/30 


Cu+Cu 


62.4 


0.32 ± 0.063 


79 


9.2 ±0.1 


14.3 ±0.2 


1.18 ±0.05 


1.92 


10"°^ ± 3.9 • 


10-03 


37.0/28 








67 


7.7 ±0.1 


12.0 ± 0.2 


1.18 ±0.05 


2.29 




10-03 


32.0/26 








57 


6.0 ±0.1 


9.1 ±0.1 


1.17 ±0.05 


2.85 


10"°^ ± 5.9 • 


10-03 


32.0/23 








48 


5.1 ±0.1 


8.1 ±0.1 


1.19 ±0.05 


3.66 


10"°^ ±8.0 • 


10-03 


29.2/21 








92 


9.1 ±0.04 


10.3 ±0.1 


1.04 ± 0.02 


4.31 


10"°^ ±9.8- 


10-04 


7.45/24 


Cu+Cu 


22.5 


0.30 ± 0.064 


58 


4.9 ± 0.02 


5.8 ±0.04 


1.06 ± 0.02 


1.11 


10"°^ ± 2.9 • 


10-03 


71.1/17 
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FIG. 3: The mean from the NBD fit as a function of TVpart 
for Au+Au (upper) and Cu+Cu (lower) collisions over the 
range 0.2 < pr < 2.0 GeV/c. The mean shown is within the 
PHENIX central arm spectrometer acceptance. The error 
bars represent the standard deviation of the distribution. 



as the fractional deviation from a scaled variance of 1.0 
of a Poisson distribution as follows: 



^ch,dyn 



1) 



(t^Gauss — 1) 



(^ch.raw 1) /gco (^ch,raw 1); 

(11) 

where Wch,dyn represents the estimate of the remaining 
dynamical multiplicity fluctuations and uJch,ia.w repre- 
sents the uncorrected multiplicity fluctuations. Since the 
correction, /geo, is calculated as a ratio of the two run- 
ning conditions of the simulation, most multiplicity fluc- 
tuations intrinsic to the model should be canceled. The 
correction always reduces the magnitude of the measured 
scaled variance. Note that the value of fgco is mathemat- 



ically identical when applied to the inverse of fcNBD : 



^NBD.dyn " /f 



gco^NBD- 



(12) 



The resulting geometrical correction factors for each 
species are constant as a function of centrality, there- 
fore a single correction factor is calculated for each trans- 
verse momentum range by fitting the correction factors 
as a function of A^part to a constant. This behavior 
is expected since centrality bins are defined to be con- 
stant percentages of the total geometric cross section. 
The correction factors for each transverse momentum 
range for a given collision species are consistent with 
each other. The standard deviation of the individual 
geometrical correction factors from the linear fits as a 
function of A^part are included in the systematic error 
of the correction factor estimation and propagated into 
the total systematic error for each point in Wch.dyn and 
^NBD.dyn- For 0.2 < pT < 2.0 GcV/c, the geometrical 
correction factors, /geo, and systematic errors from the 
fit are 0.37 ± 0.027 for 200 GeV Au-hAu, 0.33 ± 0.031 
for 62.4 GeV Au-hAu, 0.40 ±0.047 for 200 GeV Cu-hCu, 
0.32 ± 0.063 for 62.4 GeV Cu-hCu, and 0.30 ± 0.064 for 
22.5 GeV Cu-I-Cu. The extraction of the geometrical cor- 
rection factors are inherently model-dependent and are 
also dependent on the accuracy with which the central- 
ity detectors are modelled. The effect of the latter de- 
pendence has been studied by also calculating the correc- 
tion factors using constant but non-overlapping impact 
parameter distributions for each centrality bin and com- 
paring them to the correction factors using the Gaussian 
impact parameter distributions. For all px ranges, an 
additional fraction of the value of iu;ch,dyn or k^^^) dyn 
has been included in the final systematic errors for these 
quantities. The magnitude of this systematic error is 8% 
for 200 GeV Au-hAu, 8% for 62.4 GeV Au-f Au, 11% for 
200 GeV Cu-HCu, 17% for 62.4 GeV Cu-t-Cu, and 25% 
for 22.5 GeV Cu-t-Cu. A sample comparison of the scaled 
variance before and after the application of the geomet- 
rical correction factor is shown for the 200 GeV Au+Au 
dataset in Fig. [H 



IV. RESULTS 

The scaled variance as a function of the number of 
participating nucleons, A'part, over the px range 0.2 < 
Pt < 2.0 GeV/c is shown in Fig. O For all centrali- 
ties, the scaled variance values consistently lie above the 
Poisson distribution value of 1.0. In all collision systems, 
the minimum scaled variance occurs in the most central 
collisions and then begins to increase as the centrality 
decreases. In 200 GeV Au+Au collisions, this increase is 
only observed for A^part > 200. For A^part < 200 Wch.dyn 
suggests a slight decrease but is consistent with a con- 
stant value. In 62.4 GeV Au+Au collisions, the increase 
in Wch.dyn with decreasing centrality is observed only over 
the range A'part > HO. The source of the qualitative 
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FIG. 4: Fluctuations expressed as the scaled variance as a 
function of centrality for 200 GeV Au+Au collisions in the 
range 0.2 < pr < 2.0 GeV/c. Shown are the uncorrected 
fluctuations, f^ch.raw, along with fluctuations after correcting 
for the estimated contribution from geometry fluctuations us- 
ing Equation [TT] tJch.dyn. 



differences between the 200 and 62.4 GeV Au+Au col- 
lisions is not known, although some of the differences 
could be explained by the increased contribution from 
hard scattering processes at 200 GeV compared to 62.4 
GeV. Studies performed by varying the centrality selec- 
tion cuts establish that the differences are not due to 
the differences in the centrality selection algorithm. A 
similar centrality-dependent trend of the scaled variance 
has also been observed at the SPS in low energy Pb-|-Pb 
collisions at -^snn^IT.S GeV and at forward rapidities 
(1.1< j/c.m. <2.6), measured by experiment NA49 [2^ . 
where the hard scattering contribution is expected to 
be small. The Cu-|-Cu data exhibit a weaker decrease 
in the scaled variance for more central collisions. The 
62.4 GeV Cu-I-Cu scaled variance values are consistently 
above those from the 200 GeV Cu-|-Cu dataset, but the 
two are consistent within the systematic errors for all 
centralities. 

The scaled variance has been studied as a function of 
the pt range over which the multiplicity distributions 
are measured in order to determine if any significant pt- 
dependent dynamical fluctuations arc present. Results 
for several px ranges from 0.2 < px < 2.0 GeV/c down 
to 0.2 < Pt < 0.5 GeV/c are shown in Figs. [MZl In the 
absence of pT-dependent dynamical fluctuations, restrict- 
ing the Pt range should reduce the scaled variance in the 
same manner as for a fractional acceptance. Similar to 
Equation ITOl) : 



= 1 + /pt(t^rcf - 1), 



(13) 



where cup^ represents the fluctuations in the pT range of 



interest, oj^of represents the fluctuations in the reference 
Pt range, and /p.^ = UpT^/^rci is the ratio of the mean 
multiplicity in the two ranges. Also shown are curves rep- 
resenting the expected scaling of the fluctuations using 
the range 0.2 < pT < 2.0 GeV/c as the reference range. 
The shaded regions reflect the systematic errors of the 
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FIG. 5: Fluctuations expressed as the scaled variance as a 
function of A'^part for Au-|-Au (upper) and Cu+Cu (lower) 
collisions for 0.2 < pr < 2.0 GeV/c. The estimated contribu- 
tion from geometry fluctuations has been removed. Results 
from the superposition model are overlayed with the shaded 
regions representing a one standard deviation range of the 
prediction for the fluctuation magnitude derived from p+p 
collision data. Also shown (upper) is the estimated contribu- 
tion from non-correlated particle emission with the Poisson 
distribution of the scaled variance of 1.0 with the addition of 
elliptic flow in 200 GeV Au+Au collisions. 
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reference range. For all px ranges, the scaled fluctuation 
curves are consistent with the data, indicating that no 
significant pT-dependence is observed, although the data 
in the range 0.2 < px < 0.5 GeV are consistently above 
the scaled reference curves. The p^-dependence can also 
be examined in more directly with the parameter /cnbd 
from the NBD fits. Substitution of the scaled variance in 
Equation ([5]) into Equation shows that /cnbd should 
be independent of pr in the absence of pr-dependent dy- 
namical fluctuations. As shown in Figs. [5][ni there is no 
significant py-dependence of the observed values of Anbd- 
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FIG. 6: Scaled variance for 200 (upper) and 62.4 (lower) GeV 
Au+Au collisions plotted as a function of A'^part for several pr 
ranges. The lines represent the data for the reference range 
0.2 < pt < 2.0 scaled down using the mean multiplicity in 
each successive pt range. The shaded areas represent the 
systematic errors from the reference range. 
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FIG. 7: Scaled variance for 200 (upper) and 62.4 (lower) GeV 
Cu+Cu collisions plotted as a function of A^part for several pt 
ranges. The lines represent the data for the reference range 
0.2 < Pt < 2.0 scaled down using the mean multiplicity in 
each successive pr range. The shaded areas represent the 
systematic errors from the reference range. 



The scaled variance as a function of the charge sign 
of the charged hadrons is shown in Fig. [10] for 200 GeV 
Au+Au collisions in the pT range 0.2 < px < 2.0 GeV/c 
in order to investigate any Coulomb-based contributions 
to the fluctuations. In the absence of additional dynamic 
fiuctuations, the scaled variance for positively or nega- 
tively charged hadrons should be reduced from the inclu- 
sive charged hadron value by 

w+- = l + /+-(t^ch-l), (14) 
where are the fluctuations for positive or negative 
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FIG. 8: The inverse of the parameter ^nbd from the Negative 
Binomial Distribution fits for 200 (upper) and 62.4 (lower) 
GeV Au+Au collisions. The fluctuations are plotted as a 
function of A^part for several pr ranges. 



FIG. 9: The inverse of the parameter /cnbd from the Negative 
Binomial Distribution fits for 200 (upper) and 62.4 (lower) 
GeV Cu+Cu collisions. The fluctuations are plotted as a 
function of Apart for several px ranges. 



particles, Wch are the fluctuations for inclusive charged 

hadrons, and — fi^ / Hch is the ratio of the mean 

multiplicities. The scaled variance from the positive 
and negative hadrons are consistent with each other and 
consistent with the expected reduction of the inclusive 
charged hadron fluctuations. 

An additional non-dynamic contribution to multiplic- 
ity fluctuations arises from the presence of elliptic flow. 
This contribution has been estimated using a simple 
Monte Carlo simulation. In this simulation, a random 
reaction plane angle is assigned to each event. The mul- 
tiplicity distribution due to the elliptic flow component 



is given by the following: 

dN/d(j) = C[l + 2 V2 cos(2A0)], (15) 

where C is a normalization factor, V2 is the measured 
magnitude of the elliptic flow, and A(/) is the difference 
between the particle emission angle and the reaction 
plane angle. For each event, this multiplicity distribu- 
tion function is integrated over the PHENIX azimuthal 
acceptance and the resulting scaled variance from one 
million events is calculated. The value of V2 used in the 
simulation is taken from PHENIX measurements of ellip- 
tic flow at the mean transverse momentum of the inclu- 
sive charged hadron spectra in the range 0.2 < px < 2.0 
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GeV/c within the central arm spectrometers [45|- The 
estimated contribution from eUiptic flow to the observed 
scaled variance is less than 0.038% for all centralities and 
is shown for 200 GeV Au+Au collisions in Fig. [SJ 



nearly nullified after applying the geometry corrections 
described previously, so the resulting fluctuations are in- 
dependent of centrality as well as collision species. 
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FIG. 10: The scaled variance as a function of A^part for 200 
GeV Au+Au collisions in the range 0.2 < pr < 2.0 GeV/c. 
Shown are measurements for inclusive charged particles, pos- 
itive particles, and negative particles. The line represents the 
inclusive data scaled down in acceptance by 50% with the 
shaded area representing the systematic error. 



V. DISCUSSION 

A. Comparisons to a participant superposition 
model 

It is informative to compare fluctuations in relativistic 
heavy ion collisions to what can be expected from the 
superposition of individual participant nucleon-nucleon 
collisions. For this purpose, PHENIX data will be com- 
pared to a participant superposition, or wounded nu- 
cleon, model [1^ based upon data from elementary col- 
lisions. In the participant superposition model, the total 
multiplicity fluctuations can be expressed in terms of the 
scaled variance [35| . 

= ^1/ + MWN WNp„t I (16) 

where tOi, are the fluctuations from each individual 
source, e.g. from each elementary collision, c^Npart 
the fluctuations of the number of sources, and ^wn is 
the mean multiplicity per wounded nucleon. The second 
term includes non-dynamic contributions from geometry 
fluctuations due to the width of the centrality bin along 
with additional fluctuations in the number of participants 
for a fixed impact parameter. Ideally, the second term is 
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FIG. 11: Fluctuations expressed as the scaled variance as a 
function of Apart for Au-|- Au (upper) and Cu+Cu (lower) col- 
lisions for 0.2 < pt < 2.0 GeV/c. The estimated contribution 
from geometry fluctuations has been removed. Results from 
the HIJING model with jets turned on and jets turned off are 
overlayed with the shaded regions representing the systematic 
error for each curve. 

Baseline comparisons are facilitated by PHENIX mea- 
surements of charged particle multiplicity fluctuations in 
minimum bias 200 GeV p-t-p collisions. The p-l-p data 
and the NBD distribution to the multiplicity distribu- 
tion are shown in Fig. [21 The NBD fit yields ^ch = 0.32 
± 0.003, Wch = 1-17 ± 0.01, and fewBD = 1-88 ± 0.01. 



HIJING 200 GeV Cu+Cu, jets on 
HIJING 200 GeV Cu+Cu, jets off 
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These results are in agreement within errors with previ- 
ous measurements in the same pseudorapidity range of 
fcNBD = 1.9 ± 0.2 ± 0.2 by the UA5 Collaboration [20| in 
collision of protons and antiprotons at 200 GeV. Compar- 
isons of the participant superposition model to the 22.5 
GeV Cu+Cu data can be made to multiplicity fluctua- 
tions measured in 20 GeV p-l-p collisions by the NA22 
Collaboration over the same pseudorapidity range 
as the PHENIX Cu-|-Cu measurement. After scaling the 
NA22 scaled variance to the PHENIX azimuthal accep- 
tance, the participant superposition model scaled vari- 
ance is expected to be constant as a function of centrality 
with a value of 1.08±0.04. Lacking multiplicity distribu- 
tion data from elementary collisions at 62.4 GeV within 
the PHENIX pseudorapidity acceptance, it is assumed 
that as a function of collision energy, the scaled variance 
in the PHENIX pseudorapidity acceptance scales in the 
same manner as in an acceptance of 47r, which can be 
parametrized from existing p-|-p and p+p data as follows 

Mch«-4.2 + 4.69(^ri5^ (17) 

Given the mean charged particle multiplicity, the scaled 
variance in p-l-p and p+p can be parametrized as follows 

c^ch-0.35 ^^^'^~^^' . (18) 

Mch 

Scaling this parametrization to match the values of cjch 
at 200 GeV and 22.5 GeV, the estimated value of Wch at 
62.4 GeV is 1.15 ±0.02. 

Comparisons of the data to the predictions of the par- 
ticipant superposition model are shown in Fig. [5] for 
Au-|-Au and Cu+Cu collisions. The shaded regions about 
the participant superposition model lines represent the 
systematic error of the estimates described above. All of 
the data points are consistent with or below the partic- 
ipant superposition model estimate. This suggests that 
the data do not show any indications of the presence of 
a critical point, where the fluctuations are expected to 
be much larger than the participant superposition model 
expectation. 



B. Comparisons to the HIJING model 

Shown in Fig. [11] are the scaled variance curves from 
HIJING simulations into the PHENIX acceptance. The 
HIJING simulations are performed with a fixed impact 
parameter corresponding to the mean of the impact pa- 
rameter distribution for each bin as determined by the 
Glauber model in order to minimize the geometry fluc- 
tuation component of the result. The mean and variance 
of the resulting multiplicity distributions from HIJING 
are extracted from fits to Negative Binomial Distribu- 
tions. The HIJING simulation multiplicity fluctuations 



with the jet production parameter turned on are consis- 
tently above the data and increase continuously through 
the most peripheral collisions. This behavior is not con- 
sistent with the data, where the fluctuations do not in- 
crease in the most peripheral collisions. Although HI- 
JING reproduces the total charged particle multiplicity 
well, it consistently overpredicts the amount of fluctua- 
tions in multiplicity. When the jet production parameter 
in HIJING is turned off, the scaled variance as a func- 
tion of centrality is independent of collision energy, illus- 
trating that jet production accounts for the energy de- 
pendence of the HIJING results. Note that the HIJING 
results with jet production turned off are in better agree- 
ment with the data for all collision energies. Together 
with the observation that the multiplicity fluctuations 
demonstrate no significant pT-dependence, this may be 
an indication that correlated emission of particles from 
jet production do not significantly contribute to the mul- 
tiplicity fluctuations observed in the data. 

C. Comparisons to the clan model 

The clan model [i^ has been developed to interpret 
the fact that Negative Binomial Distributions describe 
charged hadron multiplicity distributions in elementary 
and heavy ion collisions. In this model, hadron produc- 
tion is modeled as independent emission of a number 
of hadronic clusters, Nc, each with a mean number of 
hadrons, rtc- The independent emission is described by a 
Poisson distribution with an average cluster, or clan, mul- 
tiplicity of Nc- After the clusters are emitted, they frag- 
ment into the flnal state hadrons. The measured value 
of the mean multiplicity, fich, is related to the cluster 
multiplicities by fich = NctTc- In this model, the cluster 
multiplicity parameters can be simply related to the NBD 
parameters of the measured multiplicity distribution as 
follows: 

Nc = A:nbd log{l + Mch/fcNBo) (19) 

and 

nc = {^J■ch/kfiBI))/log{l + Aich/fcNBo). (20) 

The results from the NBD fits to the data are plot- 
ted in Fig. [m for all collision species. Also shown are 
data from elementary and heavy ion collisions at various 
collision energies. The individual data points from all 
but the PHENIX data are taken from multiplicity distri- 
butions measured over varying ranges of pseudorapidity, 
while the PHENIX data are taken as a function of cen- 
trality. The characteristics of all of the heavy ion data 
sets are the same. The value of fic varies little within 
the range 1.0-1.1. The heavy ion data universally ex- 
hibit only weak clustering characteristics as interpreted 
by the clan model. There is also no significant variation 
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seen with collision energy. However, fic is consistently 
significantly higher in elementary collisions. In elemen- 
tary collisions, it is less probable to produce events with a 
high multiplicity, which can reveal rare sources of clusters 
such as jet production or multiple parton interactions. 
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FIG. 12: The correlation of the clan model parameters fic 
and Nc for all of the collision species measured as a function 
of centrality. Also shown are results from pseudorapidity- 
dependent studies from elementary collisions (UA5 [1711, E MC 
[ii], and NA22 ^) and heavy ion collisions (E802 ^ and 
NA35 [H). 

A feature that is especially apparent in the Au+Au 
data is the fact that the scaled variance decreases with 
increasing centrality, with the most central point lying 
below the participant superposition model expectatation. 
The clan model provides one possible explanation for 
this effect whereby there is a higher probability for con- 
tributions from cluster sources such as jet production 
in the lower multiplicity peripheral events. The cluster 
sources introduce correlations that can increase the value 
of 1/ fcNBD and hence the value of the scaled variance of 
the multiplicity distribution. Another possible explana- 
tion for this feature can be addressed with a string perco- 
lation model in heavy ion collisions \4§\. In general, per- 
colation theory considers the formation of clusters within 
a random spatial distribution of individual objects that 
are allowed to overlap with each other. The clusters are 
formed by the geometrical connection of one or more of 
the individual objects. This can be applied to estimate 
multiplicity fluctuations in heavy ion collisions whereby 
the objects are the circular cross sections of strings in 
the transverse plane [3| and the strings form clusters of 
overlapping strings that then each emit a number of par- 
ticles related to the number of strings in each cluster. As 
the centrality increases, the number of individual clus- 
ters decreases along with the variance of the number of 



strings per cluster, which can result in a decrease in the 
magnitude of the resulting multiplicity fluctuations. The 
prediction of the scaled variance from the string percola- 
tion model for 200 GeV Au-|-Au collisions scaled down to 
the PHENIX acceptance in azimuth and pseudorapidity 
[dij is shown in Fig. [TT] Although percolation describes 
the trend observed at the four highest centralities very 
well, the scaled variance from the model continues to in- 
crease well above the data as centrality decreases. The 
implementation of the HIJING model contains merging 
of strings that are in close spatial proximity, so percola- 
tion can explain the trends in the scaled variance from 
HIJING. 



VI. SUMMARY 



PHENIX has completed a survey of multiplicity fluc- 
tuations of charged hadrons in Au-|-Au and Cu-|-Cu col- 
lisions at three collision energies. The motivation for the 
analysis is to search for signs of a phase transition or the 
presence of the predicted critical point on the QCD phase 
diagram by looking for increased multiplicity fluctuations 
as a function of system energy and system volume. After 
correcting for non-dynamical fluctuations due to fluctu- 
ations of the collision geometry within a centrality bin, 
the multiplicity fluctuations in 200 GeV and 62.4 GeV 
Au-|-Au collisions are consistent with or below the ex- 
pectation from the superposition model of participant 
nucleons. The multiplicity fluctuations decrease as the 
collision centrality increases, dropping below the partic- 
ipant superposition model expectation for the most cen- 
tral Au-|-Au collisions. Fluctuations from Cu-f Cu colli- 
sions exhibit a weaker centrality-dependcncc that also is 
consistent with or below the expectation from the partic- 
ipant superposition model. The absence of large dynami- 
cal fluctuations in excess of the participant superposition 
model expectation indicate that there is no evidence of 
critical behavior related to the compressibility observable 
in this dataset. There is also no significant evidence of 
dynamical fluctuations that are dependent on the trans- 
verse momentum or the charge of the particles measured. 
As interpreted by the clan model, the observed fluctu- 
ations demonstrate only weak clustering characteristics 
for all of the heavy ion collision systems discussed here. 
The decreasing scaled variance with increasing central- 
ity may be explained by percolation phenomena, however 
this fails to explain the most peripheral Au+Au data. Al- 
though this analysis does not observe evidence of critical 
behavior, it does not rule out the existence of a QCD crit- 
ical point. Further measurements will be possible during 
the upcoming low energy scan program at RHIC allowing 
for a more comprehensive search for critical behavior. 
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